Fast computation of the Kohn-Sham susceptibility of large systems 
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For hybrid systems, such as molecules grafted onto solid 
surfaces, the calculation of linear response in time dependent 
density functional theory is slowed down by the need to cal- 
culate, in ~ TV 4 operations, the susceptibility of N non inter- 
acting Kohn-Sham reference electrons. We observe that this 
susceptibility can be calculated N times faster, with a preci- 
sion of ~ (Aoj/e) 2 , where Alu is the width of the frequency 
intervals and e a displacement off the real frequency axis. 
By itself or in combination with previously known procedures 
of accelerating such calculations, our procedure should sig- 
nificantly facilitate the calculation of TDDFT response and 
optical spectra of hybrid systems. 

Motivation. Time dependent density functional theory 
(TDDTFT) has proved to be very successful in determin- 
ing optical properties of molecules and clusters, for a re- 
view see [1]. Therefore it would be interesting, to study 
within TDFFT linear response, the simplest version of 
this theory, also hybrid materials, such as dyes bound to 
semiconductor or metal surfaces. Such materials exhibit 
interesting phenomena, for example charge transfer upon 
photo excitation in the Gratzel cell [2] or increased flu- 
orescence for organic molecules chemisorbed onto noble 
metal clusters [3]. 

Density functional algorithms based on localized or- 
bitals, such as the siesta package [4] have been used to 
study the ground state of systems as large as DNA chains 
[5], but the computation of optical spectra and proper- 
ties of excited states of hybrid systems in such schemes 
appears still to be very difficult. 

The present note is motivated by recent progress in the 
calculation of optical spectra in TDDFT linear response 
[6] using localized orbitals and the observation that even 
this recent advance is not sufficient to deal with hybrid 
systems. To model such systems one must deal with a 
slab of a few layers of solid onto which a molecule may 
be adsorbed, without the benefit of any symmetries, and 
with N, the number of orbitals, at least of the order of 
a few thousand. In the local orbital approach of [6] the 
needed CPU time grows rather steeply, as ~ iV 4 , with 
N and paradoxically, most of the CPU time is actually 
spent on the charge susceptibility of the non interacting 
Kohn-Sham reference fermions. 

Here we focus on the Kohn-Sham charge susceptibil- 
ity and describe a simple way to calculate it N times 
faster than by straightforward computation, in ~ TV 3 
rather than ~ N A operations, with a relative precision 
of ~ (Auj/e) 2 , where Auj is the spacing of frequencies, 



and we also study the performance of this technique by 
applying it to a random hamiltonian. 

The equations of TDDFT linear response [8] are based 
on the definition of the time dependent Kohn-Sham po- 
tential Vj(s( r ,t) of non interacting reference electrons: 

V KS {r,t) = V ext (r,t) + [ dr'p^ + V xc (n(r,t)) (1) 
J \r — r\ 

where V ex t is the external potential and V xc the exchange 
correlation potential which reduces to an ordinary func- 
tion of the density n(r,t) in the local and adiabatic 
approximation. Differentiating both sides of this equa- 
tion with respect to the charge density n(r,t) one finds 
a Dyson like equation for the interacting susceptibility 
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with a local exchange kernel f xc = ^jr*" 1 . The poles of 
Xint as a function of frequency provide information on the 
excited states. The preceding equation requires as input 
the susceptibility \ of independent Kohn-Sham reference 
fermions 
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where {4>i(r),Ei}, i = 1..N is a set of N Kohn-Sham 
eigenfunctions and energies. From eq (3) we see that x 
is a sum of ~ ./V 2 terms that must be computed for a 
set of N u frequencies and with r, r' each ranging over 
V/AV cells in space, with V the volume of the system 
and AV" an elementary grid volume. Because the num- 
ber of orbitals N grows linearly with the volume V, the 
calculation of Kohn-Sham susceptibility \ then requires 
~ N i N 0J operations, more than the remaining operations 
needed to calculate the interacting susceptibility Xint , the 
number of which grows with N only as ~ N^N^. 

Accelerated Computation. In LCAO schemes such as, 
for example, the siesta approach [4], one develops the 
fermion operators in terms of a set of orbitals {ip a {f)} 
of finite range. This implies a corresponding expansion 
for the Kohn-Sham density response function x an d one 
must now find 

x abcd {t) - E r a E ^r c F ^ .. n(E) ~Z {F l, - (4) 
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Here the tp^ are eigenvectors of the Kohn-Sham hamil- 
tonian in the orbital basis and the indices of Xabcd each 
range over N orbitals (which we take to be real, so we 
drop the complex conjugation). There appear to be far 
too many independent quantities Xabcd, but for orbitals 
of finite range locality limits the distance between, re- 
spectively, the orbitals (a, c) and (b,d) to a few atomic 
distances which is also quite natural as these pairs of or- 
bitals correspond to coincident points in the continuum 
limit. Therefore, in a basis of localized orbitals, only a 
subset of ~ N 2 susceptibilities Xabcd contributes in the 
linear response equations, see [6] for a more detailed dis- 
cussion. 

The Kohn-Sham susceptibility of free fermions is just 
a product (or convolution in frequency) of propagators, 

Xabcd(t) = -iG ab (t)G cd (-t) (5) 

so its calculation is essentially finished once we know the 
propagator of the free Kohn-Sham fermions. The deter- 
mination of G a h(w) requires of the order of ~ TV 3 opera- 
tions per frequency (even if the cigenfunctions are already 
available) or ~ AA W operations for a given set of indices 
{a, 6, c, d} at all frequencies, and these operations totally 
dominate the low ~ A w log A w calculational cost of a fast 
convolution. 

For this reason, and also because convolutions of Green 
functions require great care to avoid finite size effects, we 
give here a way of calculating Xabcd(u) that is compa- 
rable in overall speed, but more accurate than a con- 
volution. We decompose the susceptibility into pos- 
itive and negative frequency components, Xabcd{u) — 

Xabcdi") + Xcdab(-U) with 

E<Q<F v ' 

where pf 6 = Vf and rewrite Xabcdi^) as 

X- bcd (w) = -£ P E ab G + cd (-« + E) (7) 

_E<0 

where Gf d (u>) denotes the positive/negative frequency 
part of the propagator. We also discretize the frequency 
axis to a finite mesh of points which we choose to 
be uniform, for simplicity, and we interpolate G^(w) lin- 
early between the mesh points. Because G^ d (ui) varies 
on a frequency scale of e this interpolation introduces a 
relative error in Xabcd °f tnc or dcr of (Au/e) 2 where Au 
is the frequency spacing. The sum in cq (7) can be rep- 
resented as a convolution, Xabcd = ^ m< ^ab ® ^td(~ tJ ) 
that can be calculated much faster, but the result is less 
accurate and only useful in conjunction with an acceler- 
ated calculation of G^ b itself. 

After reconstituting Xabcd from Xabcd we nave spent 
~ NNu operations on Xabcd for a given set of indices 



{a, 6, c, d}. As mentioned before, only ~ N 2 combina- 
tions of such indices enter into the equation of the in- 
teracting susceptibility, and so we have indeed reduced 
the computational cost of the Kohn-Sham susceptibility 
from ~ N i N lx] to ~ A 3 A W operations. 

For a molecule bound to a slab of surface layers and 
similar systems that involve thousands of orbitals, this 
speed up by a factor A could turn out to be essential 
for the feasibility of the computation. 

Test of the procedure. To illustrate the quality of our 
approximate procedure as a function of the number of 
orbitals N , the number of frequency points A w and the 
smoothing parameter s, we choose a symmetric hamilto- 
nian h at random with a probability: 

prob(h) ~ exp -N £ h 2 k (9) 

i,k=l..N 

from the so called " orthogonal ensemble" of random ma- 
trices. For A — > oo, the density of eigenvalues of this 
matrix ensemble tends to a semicircular distribution in 
the interval (—1,1), see [9]. We therefore choose units in 
which the band ranges between (—1,1), set Ef = and 
calculate Xabcd (i) exactly and (ii) approximately in, re- 
spectively ~ N 2 Nu and ~ NN U operations and compare 
our results. 

With a reasonable choice of the values of e and the 
number of mesh points N u , the exact and approximate 
values of \ coincide to within the linewidth of the plot, 
as seen in figure (1). Figure (2) confirms that the error 
scales like (Aui/e) 2 with the mesh spacing. For N — 
2000, the approximate calculation is about a thousand 
times faster than the exact one. 

Conclusion. In summary, we have described an accel- 
erated procedure for computing the charge response of 
Kohn-Sham reference fermions in order to facilitate the 
calculation of TDDFT linear response for systems that 
are intrinsically of large size, such as molecules bound 
to solid surfaces. Our speedup from N^N^ to A 3 7V W 
operations for calculating the subset of elements of Xabcd 
that contribute in linear response comes at the price of a 
finite precision of the order of <~ (^-) where Auj is the 
frequency spacing. 

A reduction from A 4 iV w to A^A^ operations in 
TDDFT linear response was achieved before by limiting 
the density response to a finite manifold of fit functions 
[7]. Because our procedure is very simple it might be 
combined with such or similar methods to further reduce 
the growth of operations with A. 

We expect our algorithm to facilitate the study, by 
TDDFT linear response, of optical properties of hybrid 
systems such as molecules docked on surfaces, systems 
that are difficult to treat by existing computational meth- 
ods and which exhibit many interesting features that de- 
serve further study. 
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Figure Captions: 

Figure (1) : Rexabcd, ^Xabcd as a function of fre- 
quency for N = 2000 orbitals, e = 0.03, = 256 
frequencies and for a random choice of {a,b,c,d}. The 
difference between the exact and the approximate val- 
ues of Xabcd is smaller than the linewidth of the figure, 
therefore one can only see two distinct curves. 

Figure (2) The relative error in Xabcd for N = 2000 
orbitals and e = 0.03 as a function of the number of 
frequency points N u . The relative error is defined as the 
ratio of the root mean squares of error and signal and, as 
expected, it scales as ~ N~ 2 or ~ (Aw) 2 . 
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